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We present thorough analysis of the nonhnear tunnehng of Bose-Einstein condensates in static 
and accelerating two-dimensional lattices within the framework of the mean-field approximation. 
We deal with nonseparable lattices considering different initial atomic distributions in the highly 
symmetric states. For analytical description of the condensate before instabilities are developed, 
we derive several few-mode models, analyzing both essentially nonlinear and quasi-linear regimes 
of tunneling. By direct numerical simulations, we show that two-mode models provide accurate 
description of the tunneling when either initially two states are populated or tunneling occurs 
between two stable states. Otherwise a two-mode model may give only useful qualitative hints for 
understanding tunneling but does not reproduce many features of the phenomenon. This reflects 
crucial role of the instabilities developed due to two-body interactions resulting in non-negligible 
population of the higher bands. This effect becomes even more pronounced in the case of accelerating 
lattices. In the latter case we show that the direction of the acceleration is a relevant physical 
parameter which affects the tunneling by changing the atomic rates at different symmetric states 
and by changing the numbers of bands involved in the atomic transfer. 
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PACS numbers: 



I. INTRODUCTION 



Tunneling in weakly nonlinear systems with periodi- 
cally varying coefficients is intrinsically related to insta- 
bilities and, therefore, from the physical point of view 
the phenomenon differs fundamentally from its linear 
counterpart. This has been suggested in Ref. 1], where 
the two-level model for the Landau-Zener (LZ) tunnel- 
ing was introduced, and supported by the direct nu- 
merical simulations of the two-dimensional (2D) Gross- 
Pitaevskii (GP) equation Q where resonant tunneling 
between bands was stimulated by the nonlinearity rather 
than by the linear force and therefore termed as a non- 
linear tunneling. LZ tunneling has also been exten- 
sively studied in a number of experiments carried out 
with Bose-Einstein condensates (BECs) loaded in mov- 
ing optical lattices (OLs) fs', ^ , where in particular, it 
was found that one can distinguish two situations cor- 
responding to the instability regime and to the regime 
of Bloch oscillations. In Ref. ^ destructive effect of the 
nonlinearity in Wannier-Stark problem was investigated 
numerically. Very recently several experimental results 
were reported where resonant tunneling in tilted lattice 
was observed and the role of asymmetry of the OL on 
the tunneling was studied Q- 

The problem of tunneling has been addressed theoret- 
ically in two different limiting cases. In the case of a 
shallow lattice the tunneling was considered within the 
framework of the phenomenological models suggested in 
Q and later on thoroughly analyzed and improved in Q . 
The model accounting spatial dependence of the realistic 
field, but still based on the two-level approach was devel- 



oped in [T] . While analytical description of the crossover 
among the above limiting cases is hardly achievable, to 
the best of authors knowledge there exist no direct nu- 
merical simulations allowing one to verify validity of each 
of the models and/or to study the transient regimes, and 
no generalization of the above theories on 2D and 3D 
lattices exists as yet, in contrast to the case of the linear 
tunne ling that is well studied in multidimensional set- 
tings [UIIH. Such a generalization is the main goal of 
the present paper. 

The work is organized as follows. In Sec.|TT]we outline 
the statement of the problem, presenting the analytical 
models with different types of the lattices. In Sec. IIIII we 
deal with various cases of LZ, inter-band, and intra-band 
tunneling of matter waves in a deep OL, reducing the 
GP equation to different two- and three-level models. 
LZ tunneling in a shallow lattice is considered in Sec. lIVI 
Numerical simulations of different types of tunneling in 
2D lattices are described in Sec. |Vl The results obtained 
are summarized in Conclusion. 



II. STATEMENT OF THE PROBLEM 

We start with the dimensionless 2D GP equation 

z9tV = -AV' + V^(r,t)V' + a|V'|V, (1) 

where time t and coordinates r = {x,y) are measured in 
units 2/ijJz and = i/Zi/mLj^, respectively, with ujz be- 
ing the linear harmonic oscillator frequency, a =sign 
where as is the scattering length of the condensate. The 
dimensionless wave function '0(r) is normalized as follows 
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N = / IV'pd^r, where A'' = iV^nasNat/iz and Nat is a 
real number of atoms. The OL potential reads 

V{r,t) = 2Vo[cos^(x - aj'^) + cos^(y - a^t^) 

+2£ cos(a; — aa;^^) cos(2/ — flyt'^)], (2) 

where the amplitude Vq is measured in terms of recoil en- 
ergy and e = ei • 62. OL is assumed to be created by the 
two pairs of the counter-propagating beams with the po- 
larization vectors ei.2, which are accelerating in the direc- 
tions orthogonal to the polarization, where a — {ax,ay) 
denotes the acceleration (for description of experimental 
realization of multidimensional lattices see e.g. p^). 

For the next consideration it is convenient to make a 
substitution 



(3) 



to introduce new independent variables R = r — at^ and 
T = t, and to rewrite Eq. ([T]) in the form 

i^T* = -Aj^* + (a • R)* + T/(R)* + fxl^'p*, (4) 
where the potential is time-independent: 

V{R) = Vo [cos(2X) + cos(2f ) + ie cos(X) cos(y ) 1 . (5) 

It turns out, however that the laboratory system of 
coordinates, where R is the radius- vector, is not the most 
convenient frame for development of the theory. This 
stems from the fact that a nonseparable lattice possesses 
an elementary cell with boundaries that are rotated with 
respect to the axes of the laboratory system. This is 
illustrated in Fig. [T] where we show the OL given by 
The simple square lattice has the period ^/2tt which is 
rotated by the angle it/ 4 with respect to the laboratory 
system. In Fig. [T] (a) besides the primitive translation 
vectors li = \/2tt{1,1) and I2 = •\/27r(— 1, 1), we depict 
in the inset the basis vectors bi = l/\/2(l, 1) and b2 = 
l/V2(-l,l)ofthe reciprocal lattice which define the first 
Brillouin zone (BZ). 

Therefore it is convenient to introduce the rotated co- 
ordinate system IT] defining X = X + Y and Y = X~Y. 
In the new coordinates the potential ^ reads as 

V(R) = 2eVo (cosX + cosy + e-icosXcosr) (6) 

and represents a simple square lattice with the primitive 
translation vectors parallel to the rotated X and Y axes. 
In the new frame the GP equation reads as 

i^T* = -2Ar* + (a' • R)* + V{K)-9 + ct|*P*. (7) 

where a' = {{gx + ay)/2, (a^ - ay)/2). 

In the linear case inter-band tunneling represents a 
possibility for a particle to pass from one energetic band 
to another. In ID lattice this may occur due to the linear 
potential (the LZ tunneling), which can be interpreted in 
terms of changing the slope of the gap edges, leading to 




FIG. 1; (Color online) (a) Contour plot of the potential given 
by Eq. ([SJ with Vo ~ —0.55 and e = 0.25. In the inset we 
display the primitive vectors bi,2 of the reciprocal lattice, as 
well as the first BZ with the high symmetry points F, X, and 
M. (b) The eight lowest BZs with positions of high-symmetry 
points (only the first quadrant is shown), (c-f) Structures of 
the E{q) for diflterent amplitudes such as in (c) Vo = —0.55 
(full gap), in (d) Vo = —0.3823 (vanishing gap), in (e) Vo = 
—0.3 (overlapping of the first and second bands), in (f) Vo = 
—0.05 (shallow lattice). 



equal " effective" energies of the initial and the final steps 
(see e.g. [13 )■ In 2D and 3D lattices a particle can tunnel 
even without applying a linear force: in particular, in the 
case of a closed gap when initial and final states satisfy 
the energy and quasi-momcntum conservation laws: 



El = E2 and qi - q2 = Q, 



(8) 



hereafter Q is an arbitrary vector of the reciprocal lattice. 

Even weak nonlinearity changes the situation dramat- 
ically. First, it can couple in a resonant way two (or 
more) states, which do not satisfy the energy and mo- 
mentum conservation laws ([8]). This occurs due to pro- 
cess of four-wave mixing, which is efficient between the 
modes possessing the same group velocity. Second, the 
nonlinearity results in distinct stability properties of the 
different states, which in particular lead to asymmetric 
tunneling, i.e. to dependence of the tunneling on the 
initial conditions [U, Q . We therefore use the term the 
nonlinear tunneling in order to highlight essentially new 
features of the phenomenon. 

In this context we recall that accounting spatial depen- 
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dence of the models describing tunneling becomes highly 
relevant even in one-dimensional case where instabilities 
and nonexponential decay rate take place (see Refs. lij 
and [1, respectively). 



III. FINITE-MODE APPROXIMATION. DEEP 
LATTICES. 

A. Basic notations. 

The theory developed in the present section is re- 
stricted neither by the cos-like potential as that of con- 
sidered in ([4]) nor by the 2D case. We, however, limit our- 
selves only to the 2D case (3D generalization is straight- 
forward but more cumbersome) and, therefore, by choos- 
ing the geometry of the square lattice, we first consider 
Eq. ([7]) with a more general periodic potential V^(R): 
V(R) = y(R -I- Ij), where Ij {j = 1, 2) are the orthonor- 
mal lattice vectors. For the analysis of particular cases, 
we, however, always refer to the specific equation ([3]), (O. 

We hereafter consider the lattice to be large enough 
in each direction, what allows us to follow the standard 
practice to consider effects independent on boundary con- 
ditions by choosing the most convenient ones, specifically 
we impose the cyclic boundary conditions. 

In the analytical part of the present paper we deal 
only with the rotated system R considering Ar = 
where V = [dx.dy)- We also use the derivative with 
respect to other variables indicating them by subindices, 

e.g. Vq = (dq^,dq^). 

In the weakly nonlinear case the associated linear spec- 
tral problem 



£(/3Q,q(R) = i^a(q)V'aq(R): 



C = -2V^ 



F(R) (9) 



plays a prominent role. Here q is the wavevector in the 
reduced BZ: qx.y € [— gs, 9s] with qs — 1/2, a — 1,2, ... 
is a band number, ipaqCR) = e*'^^WQq(R), where MQq(R) 
are periodic functions Mc(q(R) = UQ,q(R -I- Ij), are the 
Bloch functions (BFs). 

We use the internal product defined by 



{u\ A\v) = ^ j u{r)Av{r)dr, 



(10) 



where A is an operator and V is the volume of the lat- 
tice. With the linear eigenvalue problem ([9]) we associate 
physically important quantities: the group velocity of a 
mode 

Vaq = {faci \ -4iV | (/Saq) = "^qEaiq) (H) 

and the tensor of the inverse effective mass M~q with the 
entrees 



(M-i),,=4- 



(V^aq I dk I iPaiq){'Paiq \ dl \ Ifaq) 



E-Ec^M) 



dqidqj 



where k,l = X,Y (see e.g. [ll] for more details). 

In the case where the lattice potential is not small, i.e. 
Vb > 1, we seek the solution of Eq. ([7]) in the form 



(13) 



where e ^ 1 is a small parameter and -tjjj are the func- 
tions of a scaled spatial, = {xj,yj) = {e^X,e^Y), and 
temporal, tj = e^T, variables (here j = 0, 1, ...). By anal- 
ogy we define Vj = [dxj ,dy-). This approach is referred 
to as a multiple-scale expansion, which in the context of 
2D BEC loaded in an OL was described in [15]. Below we 
explore some generalization of the expansion, including 
the linear force as well as resonant wave interactions in 
the consideration. 

By substituting expansion p3)) into Eq. ([7]), passing to 
the scaled independent variables, and gathering the terms 
of the same order of e, up to e'^ one obtains {j = 1, 2, 3) 

{tdt,-C)^p,^J^„ (14) 

where J^i = 0, JF2 = — {idt^ + 4VoVi)i/'i, and 

^3 = -{idt, + 4VoVi)i^2 - («9t2 + 4V0V2 + 2V?)V'i 
-I- cr]?/'! l^V'l- 

Tunneling of atoms between two states, belonging to 
the bands ai and a2 requires equality of the group ve- 
locities of the interacting modes (see e.g. [l^): 



Vi = V2 = V 



(15) 



what will be assumed in what follows. 

Finally, to shorten notations, the BFs corresponding 
to the states coupled by tunneling will be denoted as 
(j = 1,2,3) (/Ja^qj = (respectively Ma^q^- = Uj). 



B. Band structures 

Focusing on the lowest bands, in Fig. [T] we show four 
qualitatively different configurations, which can be real- 
ized by varying the amplitude of the potential. They cor- 
respond to the existence of the full gap between the first 
and second lowest bands, indicated by shaded area in the 
panel (c); closed (or infinitesimally small) gap, indicated 
by the thin horizontal line in the panel (d) ; absence of a 
full gap in a relatively deep lattice shown in the panel (e); 
and the band structure associated with a shallow lattice 
depicted in the panel (f). 

Different aspects of instabilities and inter-band transi- 
tions corresponding to the cases (c) of a full gap and (e) of 
absence of a gap, when resonance conditions of the four 
wave mixing are not satisfied, have been considered in 
[IB] and [2], respectively, and they will not be addressed 
here. We only note that before instabilities are devel- 
oped the problem is well described within the framework 
of the two-level model. 
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C. Nonlinear inter— and intra— band tunneling. 

Let us start with the case of the resonant nonhnear 
tunnehng between two different states having the same 
energy, i.e. satisfying the conditions 

El = E and qi = q2 + q2 - qi + Q- (16) 

Both p6)) and (fT5| can always be achieved in all orthog- 
onal lattices if one considers two X points related by the 
rotation by 7r/2: then qi = {qs, 0) and q2 = (0, qs); and 
in lattices with a closed gap as it is shown in Fig. [1] (d): 
then qi = ((7b,0) and q2 = {qB,qB) and the points X 
and M are considered. The most significant difference 
between the mentioned two cases is that in the first case 
the effective masses of the both states have the same sign 
(and thus the respective Bloch states posses the same sta- 
bility), while in the second case the effective masses have 
different signs (and one of the Bloch states is necessarily 
stable, while the other one is unstable). 

We look for a solution of (fT4|) with j = 1 in the form 

iPi = [Ai(ri,ii)</Ji(ro)+A2(ri,ti)(^2(ro)]e-^^*», (17) 

where Aj[xi, ti) are slowly varying amplitudes which de- 
pend on slow variables with the fastest ones being in- 
dicated explicitly. For the next step one expands (see 
Appendix rO]) 



^ i3«(ri,ti)^„q,(ro) 



^2 = 



^ Si^^(ri,ii)(^„q,(ro) 



iEtn 



(18) 



where Baq(ri,ti) are slowly varying functions (we use a 
convention that all sums with respect to a are computed 
between a — I and a = oo, with possible exclusions 
which are explicitly indicated. Substituting ^TE\\ in (fTil) 
with j — 2 and projecting on the states ipi^2 we obtain 

Ai = Ai(ri - vti,i2) and A2 = A2(ri - vti,i2). (19) 

Projection on the rest of the states (paq yields the expres- 
sion for B (see Appendix lA 2p . 

Passing to the third order we impose the orthogo- 
nality of and "01 and taking into account that the 
conditions and can be satisfied simultane- 

ously only for high-symmetric points, e.g. the points 
F, X, and M, where the group velocity is zero, v = 
and the tensor of the inverse effective mass is diago- 
nal, =diag(M-j., M'^). Then we obtain set of the 
equations 



.OA. 



- (Xiil^i 



2X12|^2|') ^1=0, 



(20a) 



-V1M2 V1A2 - X21^^? 

(X22|^2|' + 2xi2|^l|')^2 =0, (20b) 



where we define the operator 
1 92 



ViM-^Vi = 



1 92 



Ma,x dxf Ma,y dyf ' 
The nonlinearity coefficients are given by (a, P — 1,2) 

Xal3=Xl3a=Cr{fa'Pl3\(Pa'Pl3), (21) 
Xc.f3 = X0a^ ^})- (22) 

D. Nonlinear Landau-Zener tunneling. 

The LZ tunneling is induced by the linear potential 
when there exists a full gap between bands. Being orig- 
inated by the linear term it requires conservation of the 
linear momentum. In order to observe how the nonlin- 
earity affects LZ tunneling, the linear force has to induce 
an effect of the same order as the nonlinearity. Since our 
approach is valid only for weak nonlinearity, both the 
linear force and the gap, through which tunneling is con- 
sidered, have to be sufficiently small, as well. Thus we 
require 



E2~Ei^ e^£ , qi - q2 = Q, 



(23) 



where £ < 1. Further, we assume that a' = e^a with 
|a| ~ 1. Therefore, the linear force will appear only in 
the third order of the multiple-scale expansion. 

The condition together with the constraint (fT5|) 
means that in a general situation LZ tunneling can occur 
either between the two symmetric points M, with qi = 
{<lB,q_B) and q2 — {±qB,±qB), or between the points 
X, with qi (9b, 0) and q2 = (Ojigs), belonging to 
different boundaries of a gap. Therefore, as before the 
group velocity is zero, i.e. v = 0. 

Taking into account that now the energies of the two 
states are different, we perform the multiple-scale expan- 
sion for the superposition of the states 

= Ai(ri,ti)^i(ro)e-''^i*« + A2(ri,ti)^2(ro)e-^^^*° 

(24) 

and repeat the steps of derivation of (see Ap- 

pendix |X]l). Specifically, the second order term now has 
the form (here q = qi) 



iEito 



^2=E[^^q(ri,il)<^a, 

+ i?(2^(ri,ii)^„qe-^^^*° 



(25) 



The requirement of the absence of the secular terms re- 
sults in 

OT2 Z 

- (xii|^i|' + 2xi2|A2H Ai =0, (26a) 



iviM2'^ViA2 - Ke-'^*Mi 



(X22|^2p + 2xi2|^in ^2 = 0, (26b) 
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where the coupling coefRcient k is given by 

K = a((y5i|r|(^2) = a((p2|r|(/3i) . (27) 

In the case, when the Rabi frequency defined by k sig- 
nificantly exceeds the rotational frequency 5, i.e. when 
£ <^ K, the time dependence of the coupling term can be 
neglected. 

E. Interplay between Landau-Zener tunneling and 
nonlinear tunneling. 

In the cases considered above we singled out one reso- 
nant process, which is dominant for the atomic migration 
between the bands. Meantime, one can create a situa- 
tion where there exist two competing processes. This, 
in particular, may happen in the structure, depicted in 
FiglT](d), when a hnear force is appHed allowing for res- 
onant nonlinear LZ tunneling. In this particular case. 



the state M of the lowest band (state " 1" below) through 
the linear force is coupled to the state M of the second 
band (state "2" below) and through the four- wave mix- 
ing process with the state X of the second band (state 
" 3" below) . The matching conditions now read 

El = Es =E2- e^£, (28a) 

qi = q2 + Q, qi = qs - qs + qi + Q- (28b) 

To take into account all resonant processes in this case 
a three-mode model is required. Therefore, we seek a 
solution of Eq. (fT4|) in the form 

V'l = [Ai(ri,ii)(pi(ro) + A3(ri,ti)(^3(ro)]e~*^^*'' 

+A2(ri,ti)(^2(ro)e-^^^*°. (29) 

The algebra similar to one described above, results in 
the following set of equations for slowly varying ampli- 
tudes 



i-g^ + -ViMr ViAi - «e^^*M2 

- (xiil^iP + 2X12IA2P + 2xi3|v43p) Ai - xisAiAl = 0, (30a) 
+ iviM2-iViA2 - «e-^^*Mi - (x22|^2|' + 2xi2|v4ip + 2x131.431') A2 = 0, (30b) 

+ iViMg V1A3 - (X22|^2|' + 2xi2|^iP + 2xi3l^3l') A3 - X31A3AI = 0. (30c) 

I 



We mention that resonant coupling of three states (the 
two X states and one M state) is also possible in the case 
depicted Fig[T] (d) . Detailed study of the coupled- mode 
equations for the slow envelopes can be found in [l7| . 



IV. FINITE-MODE APPROXIMATION. 
SHALLOW LATTICE. 

The LZ tunneling in a shallow accelerated lattice, 
Vo ^ 1, can be most conveniently considered in the ro- 
tated coordinate system, i.e. Eq. ([7]) with the potential 
defined by Eq. ([5]). For the sake of convenience, we nor- 
malize the wave function as 5* = ip/N^^'^, and introduce 
the nonlinearity coefficient g = aN, with the assumption 
that l^l ^ 1. Now the BFs can be approximated by lin- 
ear combinations of the plane waves (see e.g. [§]), thus 
one can use the plane waves as the basis of expansion. 

We assume that the order parameter of a BEC has 
the Fourier spectrum consisting of narrow peaks, whose 
width being much smaller than the size of the BZ. In 
this case, similar to Houston's approach in the theory of 
accelerating electrons flS^, the order parameter is taken 



in the form 

^'(R,t) = e''J(*)^^Cj(i)e*'^^'^, (31) 

where the sum consists of the plane waves satisfying the 
resonant conditions ([5]), Qj are the reciprocal lattice vec- 
tors. 

Linear LZ tunneling is a result of Bragg resonance and 
occurs when the Bloch index q(t) crosses one of the Bragg 
planes (a Bragg plane in 2D is a line passing through 
the mid-point of a reciprocal lattice vector Qj, taken 
from the F-point, and perpendicular to it). The lattice 
defined by the potential (O admits two-fold (i.e. quasi 
ID) and four-fold Bragg resonances (for more details see 
Ref. (H). 

A weak nonlinearity can be accounted in the following 
way : the lattice potential and the nonlinear term are 
considered as an effective time-dependent lattice poten- 
tial Vcs = y(R)-|-(7|5'(R, t)p, where the order parameter 
is assumed periodic in R. Thus, in this approach, nonlin- 
earity does not change the relation between LZ tunneling 
and Bragg resonances. 

The potential can be rewritten in terms of the recip- 
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rocal lattice vectors as 



;Qi,oR. I jQo,iR-l 



jgjQi,iR.^_g«wi.-i".J_,_C.c. 

(32) 

Here Qi,o = (29s,0), Qo,i = (0, 295), and Qi,±i = 
Qi_o ± Qo.i (recall that qB — 1/2). The Bragg planes 
can thus be indexed accordingly, for instance, the Bragg 
plane Bi q corresponds to Qi.o- 

Two-fold Bragg resonances. The two-fold Bragg reso- 
nance occurs at any one of the Bragg planes away from 
the 0(F)-neighborhood of the high-symmetry points, 
(i.e. the points of intersection of two or more Bragg 
planes), for instance the M-points. By keeping only the 
resonant terms we obtain: 



iRi 



(33) 



The system (|35|) is identical to that in the ID lattice 
Q. The only signature of the higher dimensionality in 
the system (|35|) is that of the dependence of the sweep 
amplitude (i.e. 4qrcsa') on the angle between the Bloch 
index of the crossing point and the acceleration. 

In the linear case, g ~ Q, the tunneling probability P 
corresponding to the system ([35]) is given by the well- 
known Landau-Zener-Majorana formula p^ : 



P = |ci((X))|^ = exp ■ 



4|qrcsa'| 



(36) 



where it is assumed that "initially" |ci(— oo)| — 1. For 
instance, at the X-point we get Px = e"'^'^ ^0 /'y°-^+°-y) . 



+g(|Cip + |C2n 



gCiC2)e-''^^ 
(34) 



where Vq is the Fourier amplitude corresponding to Q. 
Substituting the expressions ([55]) and ([M)) into equation 
([7]) we get q = —a', which cancels the linear potential, 
and the set of the equations for the amplitudes Ci and C2 . 
The latter can be simplified when we use the conservation 
of the norm, in this case |Cip -I- = 1, and the phase 



transformation C, 



e-'-^c, with 



3,9/2 



(a')^t^), where we denote qros the crossing point of the 
corresponding Bragg plane. Setting q(t) = qics — a'i (i.e. 
the resonance is set at i = 0), we obtain the set of the 
equations (j = 1, 2) 



4qr 



.a't+|(|c,f 



(35) 



Four-fold Bragg resonances. This resonance occurs at 
the intersection of two Bragg planes, for instance at the 
M-point [the other case is at the F-point, see FiglT] (f)]. 
Using the lattice potential in ((32| . the resonant part of 
the order parameter 



(37) 



^Q(t)e'('J^"'^'*)^, 

where qi = Qi,i/2, q2 = Q_i,_i/2, qs ^ Q_i^i/2, and 
q4 = Qi^_i/2 connect the F and M points, and the re- 
lations Qi^ia' — Ox and Qi._ia' ~ Oy we obtain a set of 
the equations for the Fourier amplitudes Cj of the order 
parameter. The system is simplified after the transfor- 



mation Ci 



with 



using of the norm conservation ^ 



- 2[q2 + (a')2i2] 



and 



icx = (-2aa;i - g|cip) ci + -yC2 -I- £^0(03 + C4) + 2gc2C3C4, 

ic2 = (2a2:i - g|c2p) C2 -I- -yci -f eVo(c3 -f C4) -I- 2.gciC3C4, 

ic3 = (2ayi - g|c3p) C3 + -yC4 + eVo(ci -I- C2) -f 25C1C2C4, 

ici = (-2ayt - .g|c4p) C4 -I- -yC3 -f eVo(ci + C2) + 2gciC2C3. 



(38) 
(39) 
(40) 
(41) 



In the linear case, g = 0, one can use the general for- 
mula for the tunneling probability [131, P = |cj(oo)p = 

exp{-27rX:„,^j- |[t'^\!'| } for \cj{-oo)\ = 1, where ly^ 
is the strongest sweep amplitude and A^j is the cross- 
coupling term. For instance, in the case of \ax\ > \ay\ we 



get for P = |ci(cxd)|^ 



P = exp < — 



nVl 



(42) 



with the initial condition |ci(— oo)| = 1. The same for- 
mula is true for P = |c4(oo)p with |c4(— cx))| = 1. The 
physical meaning of Eq. (|42p is clear: it determines the 
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fraction of BEC atoms which pass through the border 
of the first BZ at the M-point by continuously increas- 
ing their quasimomentum (BEC atoms leave the first BZ 
also through the side M-points, with the corresponding 
fractions given by the Fourier amplitudes C3 and C4). 

In the special case of a^; = = a using formula p6p . 
we obtain the probability of the transition in the form 



[defined after ((T])] what corresponds to 3 x 10'* and 6 x 10^ 



of real atoms. 



A. Stationary optical lattice 



Non-tunneling regime and the modulational instability 



P — exp < — TT 



2\a\ 



(43) 



Finally, the squared amplitudes |cj(t)p of the LZ sys- 
tem can be compared with the powers of the Fourier 
peaks in the full PDF (defined as the integral of the 
modulus squared of Fourier amplitude) - see Fig. [Tl] 
of Section 



V. NUMERICAL RESULTS 

The theory of tunneling developed in the preceding sec- 
tions is based on the two- (few ~) mode approximation. 
As it has been shown in 2] on an example of inter-band 
tunneling between two initially populated states, this ap- 
proach is valid until instability is developed. In general, 
when initially only one of the states is populated or when 
one assumes several populated states, that are not con- 
nected by some of the resonant conditions, the theoretical 
model requires verification. 

In this Section we address problems of tunneling start- 
ing with only one populated state in a stationary lat- 
tice and of LZ tunneling in an accelerating lattice with 
different directions of the lattice acceleration and differ- 
ent lattice depths. We impose the initial conditions in 
the form ipi^^ 0) = X]q q "^"q' where summation is 
taken up to the eight lowest bands, a = 1,...,8, and 
over the symmetric points which we denote as qr = 
(0,0), qx = (±1/2, ±1/2), qx' = (±l/2,Tl/2), and 
Qm = [(0,±1), (±1,0)]. We study evolution of the pos- 
itive quantities r^q representing populations of the re- 
spective states. Also we use the notation Yq for a point 
Y=(X, M, F) of a-th band. 

In all numerical simulations we solve the original GP 
equation (P), i.e. our results are obtained directly for 
the accelerating lattice ([2|). We use the Crank-Nicolson 
scheme with periodic boundary conditions on the grid 
with 400 X 400 points and with 20 periods of the OL in 
each direction (such a size of the lattice lattice is suffi- 
ciently large to allow one to observe volume phenomena 
that are not affected by the specific type of the boundary 
conditions). Finally, we assume inter-atomic interactions 
to be repulsive, i.e. a = 1. 

For a typical experimental setup with lO"' atoms hav- 
ing scattering length of order of 0.01 nm (achievable by 
means of Feshbach resonance) and with transversal size 
of the cloud £z ~ 0.3 /im our numerical unit of time cor- 
responds approximately to 0.25 ms. Also, in what follows 
we use the normalization constants = 10 and A = 20 



We start with examining the case where energy of the 
initially populated Mi-point is the same as that of X2- 
point which is not populated at t = 0. This is the case 
described by p^ . Now in contrast to the situation de- 
scribed in [2], now it is a priori not known which mode 
(or modes) the energy will be transferred to. Thus one 
may expect tunneling between bands, induced by modu- 
lational instability (MI) in the state Mi (since the effec- 
tive masses Mi^x and Mi^y are negative and a = 1). The 
results are shown in Fig. [21 
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FIG. 2: (Color online) Development of MI starting in the 
Mi-point. In (a) temporal behavior of the density maximum 
[at t « 1300 (what corresponds 0.325 s, i.e. to experimentally 
feasible time) the MI starts to develop] and in (b) the dynam- 
ics of the population of the Mi-state is shown (populations 
of the higher bands in the F-, X- and M-points are smaller 
than 10^^). The parameters are as follows: Vb — —0.3823, 
£ = 0.25, iV = 10. In (c), (d) density profiles and in (e), (f) 
their Fourier transformations at different times are shown. 

The most relevant result stems from the fact that even 
after the MI is developed the populations of the sym- 
metric states of the upper bands, including the states 
X2 are negligibly small. It means, that in spite of the 
resonant conditions (|16p being satisfied, tunneling does 
not occur, i.e. the expected scenario, which is described 
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above, is not realized. In order to investigate the dynam- 
ics in the Fourier space, we have computed the Fourier 
transform of the wave function at different moments of 
time see Fig. [2] (e), (f). In the process of the develop- 
ment of MI the atoms spread out along the boundary 
of the first BZ. Dispersion spreading proceeds until the 
MI is developed. Only the unstable state Mi imprints 
its signature in the symmetry of the developed structure 
[see Fig. [2] (d)] and therefore its pattern closely resem- 
bles the structures emerged from the unstable states in 
a separable lattice (c.f. with Fig. 7 in the Ref.jT^]). 



2. Block intra-band tunneling at the X2-point 

We obtained perhaps even more surprising results 
when in the lattice with the band structure depicted in 
Fig. [1] (d) we initially populated a stable state X2 that 
possesses two positive effective masses. In Figs. [3] and S] 
one can observe exchange of particles between the two 
doubly degenerate X2-points which are related by the ro- 
tation 7r/2 and referred to as X2 and below. In other 
words the system develops intra-band tunneling instead 
of the inter-band one. 




FIG. 3: Time dependence of the maximum of the density 
amplitude [(a),(c)] and populations r2x (solid line) and T2x' 
(dashed line) [(b), (d)] of the X-states of the second band in 
the case when in (a), (b) initially all particles are in the X2- 
state and in (c), (d) initially r2x = 0.99 and r2x' ~ 0.01. In 
(a) at t ^ 9000 (« 2.3 s) and in (c) at t ^ 7300 (^ 1.8 s) the 
MI starts to develop. The ratio between the time intervals 
when the particles populate X- or X'-states in (b) is ~ 1.2. 
The parameters are as follows: A ~ 0.792, Vo ~ —0.3823, 
e = 0.25, N = 10. 
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FIG. 4; (Color online) In (a)-(c) density profiles and in (d)-(f) 
their Fourier transformations at difi'erent times corresponding 
to Fig|3](a), (b) are shown. 



This picture agrees with the theory developed in 
Sec. nil CI Indeed, first we recall that the both states 
X2 are stable, while the state Mi having the same en- 
ergy is unstable. Next, we consider pp)) and take into 
account that initially the envelopes Ai^2 = ^1.2(^2) are 
independent on spatial coordinates: Vi^i 2 = 0. This 
allows us to reduce the system of PDEs to a Hamiltonian 
system with the Hamiltonian 



H, 



eff 



2xi2|Ai|2|yl2| 



-X12AIAI+X12AIAI. (44) 



Here we have taken into account the symmetry of the 
lattice resulting in xii = X22 and X12 = X12 = X21 = X21 
[see Eqs. dH]) and P^ ]. 

Since the derived dynamics preserves the "total num- 
ber of particles" Af = Ni + N2, where Afj = \Aj\'^ we in- 
troduce the population imbalance z = {\A2\^ — |^i|^)/A/', 
the relative argument (j) = 2a,Tg{AiA2), the param- 
eter A — X21/X111 as well as the renormalized time 
T = Xii-^^2- Now the dynamical system can be rewritten 
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in the form 
dz 



= A(l - z^) sin( 



dH, 



eff 



^ = 2z(l-Acos<^-2A) = -^, 
dr oz 

where the effective Hamihonian has the form 
Heff = Acos^-2A). 



(45) 
(46) 

(47) 



For the next consideration one should keep in mind 
that the dependent variables are bounded: — 1 < z < 1 
and < (/) < 2iT. From Cauchy inequality (I'^aPlv'/Jp) < 
^(|(pa|4)(|(p^|4), it follows that < A < 1, what in the 
case at hand is transformed into Xi2 < Xii- 

It is important to note here that the only nonlinearity 
contribution controls the scaling of the time r, while the 
parameter A is completely determined by the lattice it- 
self. Thus the nonlinearity sets the time scale, but not 
the type of the dynamics (see below) , which is defined by 
the lattice. 

The Hamiltonian (|47|) has two fixed points Fi = {z — 
0, (/) = 0} and F2 — {z — Q^cf) = tt}: Fi is a saddle point 
for A < 1/3 and a local minimum for A > 1/3, while 
F2 is always a local maximum (since A < 1). In Fig. [5^ 
we show the phase portrait corresponding to regime for 
A > Ac = 1/3. 
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FIG. 5; In (a) phase portrait of the system (|45}, (|46)) with 
A = 0.79 > Ac is shown. In (b) and (c) the oscillations of the 
normalized population imbalance with the initial conditions 
z(0) = 0.98, <^(0) = in (b) and (^(0) = tt in (c) correspond- 
ingly. Here solid and dashed lines correspond to r+ and r_. 

In the panels (b) and (c) of Fig[5] we present the dy- 
namics of the "trajectories", r±{t) = [1 ± z(<)]/2, that 
correspond to the populations r2x and r2X' in the full 



picture illustrated by Fig. [3l When all particles are ini- 
tially in the X2-state one obtains {^(0) = 1,0(0) — 0}, 
what corresponds to a saddle point. This fact explains 
why the initial dynamics of the population FigO (b) is 
constant. By time, however numerical errors accumu- 
late what results in motion of a system along the sep- 
aratrix which separates two adjacent regions with dif- 
ferent dynamics. In the course of evolution the system 
having reached the next saddle point jumps from one 
region to the other where the period of oscillations is 
different. This process continues and it demonstrates 
itself through the transitions in dynamics of the popu- 
lations in Fig[3] (b). In order to verify that the transi- 
tions occur due to the fact that the system stays close 
to the separatrix we impose a small shift by taking ini- 
tially r2x = 0.99 and r2X' = 0.01, which corresponds to 
{z(0) = 0.98, 0(0) = 0}. As one can see in Fig. |1 (d) the 
system immediately starts to oscillate with the fixed pe- 
riod. This can be interpreted as the oscillations around 
the fixed point Fi [c.f. Fig. [5] (b)]. By comparing the 
ratio between the periods in different regions near the 
separatrix, in Fig[5](b), with the ratio between different 
times in Figl3] (b), we found that they are 1.4 and 1.2 
respectively. 

Returning to Fig. [31 one can observe that the particle 
transfer due to the intra-band tunneling is accompanied 
by small oscillations of the density amplitude. After some 
time the MI is developed. From the first sight this dy- 
namics looks surprising because the energy interchange 
occurs between two stable states. However, by inspect- 
ing the initial state in the points of high symmetry of the 
first five bands we found that due to inherent inaccuracy 
of the numerical orthogonalization procedure initial data 
have nonzero contributions (~ 10"'*-;- 10~^) from the un- 
stable points of higher bands (for example, in the case of 
FigEKa), (b) we have also contribution from the unstable 
points rgx' « 1.3 x 10"^ and rgx « 1.2 x 10"^). These 
small contributions of unstable states are consistent with 
the estimation of the time of MI development. 



B. Accelerated deep optical lattice 

1. LZ tunneling starting from a stable state. 

Now we turn to the numerical study of tunneling in an 
accelerating OL. Along this Section we consider acceler- 
ation of a deep OL having the band structure depicted 
in Fig[l](d). We start with the case when initially only a 
state Xi is populated and the acceleration is applied in 
the x-direction, a = [ax , 0) . The results are summarized 
in FiglS] (a), (b) where the populations are obtained at 
instants when the running OL having passed an integer 
number of period returns to its initial position. In all 
our simulations we observed that the high symmetrical 
F- and M-states do not become significantly populated 
during the dynamics. 

Considering acceleration applied in the x-direction in 
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FIG. 6: Dynamics of populations of the lowest bands in X- 

5 X 10""', a-u = in 



points with applied acceleration 

5 

V2 



[(a),(b)] and with a^, = = ^ x 10 in [(c), (d)]. In 



[(a),(c)] the linear, a = 0, and in [(b), (d)] the nonlinear, 
(7 = 1, regimes with TV = 20 are shown. Different symbols 
stand for populations of X-points (the lines are guides for an 
eye). 



the linear case {a = 0) in Eq. ^ [see Fig. [S] (a)], we 
observe three different regimes. First, during the initial 
stage the atoms tunnel into the second band implying 
increase of the population of the X2-state. This regime, 
at t R:! 60 is changed to balanced distribution of atoms 
between Xi and X2 states, which lasts until t w 150. 
This stage of the evolution can be explained due to the 
large energy difference between the states X2 and Xa ^.s^g- 
In our numerics the evolution is ended up in the parti- 
cle transfer to the upper states. We inspected the eight 
lowest bands and found that the process is accompa- 
nied by developing rather complex spatial structures of 
the atomic distribution. Since the states X3^4,5_e have 
very close energies, their populations reveal almost syn- 
chronous behavior, as this is illustrated in Fig.[S](a). 

The nonlinear LZ tunneling has two additional fea- 
tures [see Fig. [6] (b)], compared to the linear one. The 
initial stage of the tunneling from Xi to X2 is ended 
at earlier times, t « 40, and during the second stage 
the dynamics reveals significant Rabi oscillations between 
the states. Such oscillations occur due to the cross- 
phase-modulation terms [see the two-state models (|20|) 
or (|26p ] and are characterized by a characteristic period 
tfiMbi ~ 100. Taking into account that now the coupling 
coefficient k ([27]) computed between the two X-states 
K « 0.011 we find that the analytical estimate for the 
same time following from the two-state approximation 
yields 1/k « 90, which is sufficiently good estimate, tak- 
ing into account that, strictly speaking, the two-state 
approximation is not valid any more. In the third stage 
one observes tunneling of the atoms to the upper bands. 

In order to understand how tunneling is affected by the 
direction of acceleration, we show in Fig[n](c), (d) the re- 



sults for flj; = ay. We obtained them in such a way that 
the total modulus of the acceleration is the same as in 
FigEl (a) , (b) . One can observe the two following changes 
in the tunneling process. First, exchange between the 
first and second bands is much stronger: during the sec- 
ond stage of evolution, 120 ^ t < 200, the second band 
acquires larger population than the first band. Second, 
after the instability is developed at t > 200 the atoms 
start to populate 7-th and 8-th bands while 3 - 6-th bands 
remain negligibly populated. 

This last phenomenon can be explained by the symme- 
try of the system. Indeed, let us consider the lowest eight 
BZs shown in Fig. [1] (b). Departing from the Xi point 
one finds the closed high-symmetry points located along 
the direction (1,1), which coincides with the direction of 
acceleration, to be the Xy^g points, while along the direc- 
tion (1,0), i.e. along the direction of acceleration shown 
in Fig. [6] (a), (b), one finds the closed points X3_4,5_6. 



2. LZ tunneling starting from an unstable point. 

Let us consider now the situation when initially all par- 
ticles are placed in the unstable point Mi. In Fig[7]we 
present linear and nonlinear dynamics of the populations 
of the lowest bands when the acceleration is applied in 
(1,0) and in (1,1) directions. In the panels (a) and (b) one 
observes that the three lowest bands participate in the 
exchange of particles, which can be explained by the fact 
that the energy differences between the states Mi, M2 
and M2, M3 are of the same order [see Fig. [T](d)]. More- 
over, the population of the third band becomes larger 
than the population of the two lowest band. 
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FIG. 7: Dynamics of populations of the lowest bands in the 
M-point subject to the acceleration (ax,ay) = (5 x 10"'', 0) 
with A'' = 20 in [(a),(b)] and to — Oy = 5 x 10"'' with 
iV = 10 in [(c), (d)]. In [(a),(c)] linear cr = and in [(b), (d)] 
nonlinear a = 1 regimes are shown. 

In the linear case the third stage of the evolution has 
not been reached for the time t < 250 (this time was 
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sufficient for transfer to upper bands in all other cases 
considered so far). Thus the effect of the nonlinearity 
is much more pronounced, than in the case of tunnel- 
ing from the stable X-point. Now the nonlinearity gives 
rise to the Rabi oscillations between the first and the 
third bands accompanied by strong attenuation of the 
population of the second band. One can interpret this 
phenomenon by means of the definition of the nonlinear 
coefficients Xq/3 [see (|211) ]. which are responsible for the 
exchange of atoms between a-th and /9-th bands, and 
from the parity of the BFs bordering a gap: BFs of the 
bands with different parity have also different parity and, 
therefore, atomic transition induced by the nonlinearity 
between first and third bands is more favorable than the 
oscillations between first and second bands. 

The results shown in Fig.[7](c),(d) indicate that change 
of the direction of the acceleration results in significant 
variations in dynamics of populations. Namely, as one 
can see in Fig 17] (c), where the direction of acceleration 
coincides with the direction of the symmetry axis (1,1), in 
the linear case there occurs a monotonic transfer of atoms 
from the first band to the upper bands, populating almost 
equally M points in the third and forth bands, with tun- 
neling to the upper bands at later times. Remarkably, 
the nonlinearity does not affect much the tunneling in 
this case [c.f. Figs. [7] (b) and [7] (d)], giving rise only to 
relatively weak oscillations of the populations. 



LZ tunneling with broken symmetry between X and M 
states. 



Figl9] is significantly larger than in the case when parti- 
cles tunnel from X- to M-state presented in FiglS] what 
reveals extremely strong effect of the instability which at 
very early times results in equal redistribution of atoms 
between the two high-symmetry points. 
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Now we focus on the situation when initially two states 
with different symmetry occupied, namely, X- and M- 
states. Specifically, in the case of the resonant tunnel- 
ing presented in Fig[l](d), we consider the two following 
configurations: i) almost all particles are in the Xi -state 
while only small number of particles are in the state Mi 
[see Fig. |S] and ii) almost all particles populate the state 
Ml while only small number of particles are in the Xi- 
state [see Fig. ^ . 

In order to verify whether there exists tunneling be- 
tween X- and M-states in we calculated the sum of the 
populations of the eight lowest bands, X]ct''x,M, in each 
symmetrical point in the linear and nonlinear regimes. 
The results are presented in Figs. [5] (e) and [S] (e). The 
thick dashed and thin dashed lines which represent pop- 
ulations in X- and M-points indicate that particles in 
X- and M-points conserve their initial population and 
thus behave independently. It means that in the case 
of broken symmetry the pure acceleration couples differ- 
ent bands but does not couple different symmetry points. 
In the presence of the nonlinearity (thick solid and thin 
solid lines) the evolution of populations reveals Rabi os- 
cillations of populations between X- and M-points which 
occur due to cross-phase-modulation. We note that this 
effect without acceleration was observed in 0| . 

The tunneling from the M- to the X-state, shown in 



FIG. 8: In (a), (b) linear, cr = 0, and in (c), (d) nonlinear, 
a = 1, dynamics of the populations of the lowest bands in X- 
and M-points with acceleration {ax,ay) = (5 x 10~*,0) with 
= 20 are presented. In the panels (a) and (c) populations 
in the X-point and in the panels (b) and (d) populations in 
the M-point are shown. Initial distribution of the particles 
is rix = 0.9 and riM = 0.1. Dynamics of the sum of the 
populations in the X- (thick solid and dashed lines) and in the 
M-points (thin solid and dashed lines) in the linear (dashed 
lines) and nonlinear (solid lines) regimes is shown in (e). 



C. LZ tunneling in a shallow lattice 

Now we turn to the case of the shallow OL correspond- 
ing to Fig[T] (f). We start with situation when all par- 
ticles populate the Xi-point. By accelerating the OL in 
the linear regime one can see in FiglTU] (a) , (b) that small 
number of particles succeeded to pass to the upper bands 
and after some time interval (in our case t « 300) the 
population of the first five bands becomes stable. The 
peculiarity of this process stems from the fact that after 
t « 200 the population of the third and fourth bands 
starts to exceed the population of the second band. This 
effect can be explained by the difference in the effective 
mass of these bands by taking into account band struc- 
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FIG. 9; The same as in Fig[8]with initial distribution of the 
particles riM = 0.9 and rix = 0.1. 



ture from FiglT] (f) the first, third and fourth bands have 
the negative effective masses while the second and fifth 
bands have the positive ones. 




0.03 



periods of oscillations this dynamics evolves into the MI. 

The dynamics in the shallow lattice case can be inter- 
preted also in terms of the Fourier powers instead of the 
Bloch band populations. Indeed, as the theory outlined 
in the Scc. lIVI predicts for the times smaller than the time 
of the instability development, the X-point corresponds 
to the two-fold quasi ID resonance, hence predicting only 
two significant peaks in the Fourier space. Namely, if 
we assume that the initial condition is broad enough in 
the real space, e.g. a Gaussian possessing width cover- 
ing many lattice periods, so that the Fourier transform 
will be a narrow peak, then one can use the LZ model 
([55]) (we notice, that while Gaussian distribution of the 
density can be substituted by any other one of the same 
spectral width, it can be created experimentally by first 
loading the condensate in a parabolic trap at low den- 
sity, or alternatively sufficiently small wavelength, and 
subsequently switching on the optical lattice with simul- 
taneous switching off the parabolic trap) . 

We observe excellent agreement in the linear tunnel- 
ing regime and a good quantitative correspondence be- 
tween the LZ models and the full PDE. For instance, in 
Fig. [11] the numerical results and theory for tunneling at 
the M-point in the linear case are compared in the case 
of four-fold resonance. One should keep in mind that 
for the initial value of the Bloch index we take into ac- 
count a time shift in the LZ model. In the nonlinear case, 
however, the theory poorly corresponds to the numerical 
results, though there is excellent qualitative correspon- 
dence easily observed in the Fourier space (for instance, 
in the number of significant Fourier peaks). 
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FIG. 10: Dynamics of projections at the X-point with applied 
acceleration a^; = 5 x 10~* in (a) the linear case, cr = and 
in (b) the same as in (a) with the presence of nonlinearity, 
(7 = 1. Here A*' = 10. Other parameters are Vb = —0.05 and 
£ = 0.25 [the band structure is shown in Fig[l](f)]. Different 
symbols denote the projections to the a-bands at the X-point 
calculated at time when the OL returns to the initial position 
(thin lines are guides for an eye). 



FIG. 11: Comparison of the LZ model with full PDE in the 
linear case. We use the integrated powers of the Fourier 
transform of the order parameter ^ (solid lines) and the 
squared absolute values of the LZ coefficients (dashed lines). 
The horizontal line indicates the theoretical prediction ac- 
cording to formula (fi^. Here Vb = 0.025, e = — 1 and 
(aa;,ay) = (0.01,0). The Gaussian initial condition for \E' 
with width of eight lattice periods are used. 



One observes completely different dynamics in the non- 
linear case with the same initial conditions as in the linear 
one. In FiglTO] (c) the population of the particles in the 
presence of nonlinearity starts to oscillate between the 
first and the second bands and only very small number 
of particles transfers into the upper bands. After several 



VI. CONCLUSION 

In the present work we have carried out extensive an- 
alytical and numerical analysis of the nonlinear tunnel- 
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ing of Bose-Einstein condensates in static and accelerat- 
ing two-dimensional lattices. We focused on nonsepara- 
ble potentials possessing different band gap structures: 
having a full gap, no gap, and satisfying the resonant 
conditions with the full gap vanishing. Wc investigated 
different initial atomic distributions in the highly sym- 
metric states, which in the nonlinear case can be either 
stable or modulationaly unstable, showing substantially 
different evolutions. We also studied the effect of the 
direction of the lattice acceleration on the nonlinear tun- 
neling. Within analytical description of the process, we 
developed several two- and three-mode models, assuming 
that only a few states are occupied by the atoms. 

We found that two-mode models accurately describe 
the tunneling before instabilities are developed, in the 
cases when either two respective states are initially pop- 
ulated or tunneling occurs between the two stable states, 
as this is the case of intra-band tunneling. The early 
stage evolution is dominated by the Rabi oscillations. 
Otherwise the two-mode models, although giving useful 
qualitative hints for understanding of tunneling and even 
estimates the characteristic times of the process, in gen- 
eral do not describe adequately all the peculiarities of 
the phenomenon. This deficiency can be explained by 
the fact that we deal with spatially distributed systems, 
which in the nonlinear case develop instabilities induced 
by the two-body interactions, which in their turn result 
in emergence of spatial coherent structures which are 
associated with population of higher bands. Moreover, 
a two-dimensional band structure possesses such pecu- 
liarities as states having tensor of the inverse reciprocal 
masses with components of different signs, what natu- 
rally diversifies dynamical scenarios of the tunneling and 
the symmetry of the developed patterns. We also found 
that in accelerating lattices the instability phenomenon 
is much more pronounced and is developed at relatively 
early stages of evolution. In the latter case we have shown 
that the direction of the lattice acceleration is a relevant 
physical parameter significantly modifying atomic pat- 
terns emerging as a result of interplay between tunneling 
and instabilities. 
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APPENDIX A: DETAILS OF MULTIPLE-SCALE 
EXPANSION 

1. Some useful relations 

In this paper we use the relations as follows 

(y'aiqi |r|¥'Q2q2) = ("aiqi l^l ^Qaqs >'5qi ~q2 ,Q i (^1) 



ii|Vq|u2) = -i{ipi\r\ip2) = -4i 



. (<y9l|V|(^2) 



E2 — El 



(uilVlua) 



(A2) 



E2 — El 

(A3) 

21 being valid for a periodic BFs bordering gap edges. 



2. Derivation of Eqs. ((20)) . 

The higher order terms in the expansion (|13p are 
searched to be orthogonal to the leading order term ipi: 
if 1,2 I ipj) = U = 2, 3). Therefore we represent '02 as 



4'2 = y^^Baq{ri,ti)ipaqe~ 



iEto 



(A4) 
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and choose the coefficients Baq{ri,ti) in the form 

Baq = -B£^'(5q-qi,Q(l — Saai) + -S^^' (Jq-q^ .Q (1 - (JcqJ, 

(A5) 

what justifies the formula 

In order to obtain (dH) we substitute (|A4p with (|A5|) 
into what yields 



^ K(qi)-i?]B«^„q, + ^ [i?„(q2)-i?]i?i2Vaq. 

= idt.Aiipi + idt,A2Lp2 + 4ViAi • Vo^pi + 4V1A2 • Vot/?2. (A6) 



By applying (931,2 | to the both sides of this equation and using (IA3|) we obtain 

dt.Aj - vViAj = for j = 1, 2. (A7) 



14 



Next we apply (i^a^ .q^ ^ \ into (|A6p and use (jA3[) and ([T9l 
This yields for j = 1 , 2 and a ^ ofj 



^0-) ^ iVol^,) 



(A8) 



The equations of the third order (PO)) are obtained from 
(|9|) by using the orthogonality conditions (<y5i,2 | -^s) = 0, 
the explicit form of ?/'2 given by and (|A8[) . and the 
definition (fl! 



3. Derivation of Eqs. pS]) . 

Now the coefficients i?i"'q(ri, ii) are as follows 

S^ci= Vq,.Q(l-'5ac.,)i3i-''\ J = 1,2 (A9) 

what justifies the formula 

In order to obtain ([19]) we substitute ([251) with (|X9l) 
into (|14p . separate terms corresponding to the main har- 



monics [cx exp(— ii?i^2io)] what yields (j = 1,2) 

-4ViAj • Vo(^j. (AlO) 

By applying (<y5i,2| to the both sides of this equation 
and using (fTT|) we compute Eqs. (|A7|) . Next we apply 
(V'aj,qi.2l to (PlOl) and use and (HH) . This yields 
for j ~ 1,2 and a ^ 



(All) 



The equations are obtained from ((51) after sepa- 
rating the two main harmonics, using the orthogonality 
conditions (1^91.2 | ^3) = 0, a particular explicit form of 
-02, dsn]), and the definition P^ . 



[1] V. V. Konotop, P.G. Kevrekidis, and M. Salerno, Phys. 

Rev. A 72, 023611 (2005). 
[2] V. A. Brazhnyi, V. V. Konotop, and V. Kuzmiak, Phys. 

Rev. Lett. 96, 150402 (2006). 
[3] O. Morsch, J. H. Muller, M. Cristiani, D. Ciampini, and 

E. Arimondo, Phys. Rev. Lett. 87, 140402 (2001); M. 

Cristiani, O. Morsch, J. H. Muller, D. Ciampini, and 

E. Arimondo Phys. Rev. A 65, 063612 (2002); C. Sias, 

A. Zenesini, H. Lignier, S. Wimberger, D. Ciampini, O. 

Morsch, and E. Arimondo, Phys. Rev. Lett. 98, 120403 

(2007). 

[4] M. Jona-Lasinio, O. Morsch, M. Cristiani, E. Arimondo 

and C. Menotti, 'cond-mat/0501572 " 
[5] S. Wimberger, R. Mannella, O. Morsch, E. Arimondo, 

A. R. Kolovsky, and A. Buchleitner, Phys. Rev. A 72, 

063610 (2005). 

[6] C. Sias, A. Zenesini, H. Lignier, S. Wimberger, D. 
Ciampini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 
98, 120403 (2007). 

[7] G. Ritt, C. Geckeler, T. Salger, G. Cennini, and M. 
Weitz, Phys. Rev. A 74, 063622 (2006); T. Salger, G. 
Ritt, C. Geckeler, S. Kling, and M. Weitz, (in prepara- 
tion). 

[8] B. Wu and Q. Niu, Phys. Rev. A 61, 023402 (2000); O. 
Zobay and B. M. Garraway, Phys. Rev. A 61, 033603 
(2000). 

[9] V. S. Shchesnovich and S. B. Cavalcanti, J. Phys. B: At. 
Mol. Opt. Phys. 39, 1997(2006). 
[10] V.S. Shchesnovich, S. B. Cavalcanti, J. M. Hickmann, 
and Yu. S. Kivshar, Phys. Rev. E 74, 056602 (2006); A. 
S. Desyatnikov, Yu. S. Kivshar, V. S. Shchesnovich, S. 



B. Cavalcanti, and J. M. Hickmann, Optics Lett. 32, 325 
(2007). 

[11] A.R. Kolovsky and H.J. Korsch, Phys. Rev. A 67, 063601 

(2003) ; D. Witthaut, F.Keck, H.J. Korsch, S. Mossmann, 
New J. Phys. 6, 41 (2004). 

[12] M. Greiner, I. Bloch, O. Mandel, T. W. Hansch, and 
T. Esslinger, Phys. Rev. Lett. 87, 160405 (2001); M. 
Greiner, O. Mandel, T. W. Hansch and L Bloch, Nature 
419, 51 (2002); H. Moritz, T. Stoferle, M. Kohl, and T. 
Esslinger, Phys. Rev. Lett. 91, 250402(2003). 

[13] J. M. Ziman The principles of the Theory of Solids (Cam- 
bridge: Cambridge University Press, 1972). 

[14] P. Schlagheck and S. Wimberger, Appl. Phys. B 86, 385 
(2007). 

[15] B. B. Baizakov, V. V. Konotop, and M. Salerno, J. Phys. 
B: At. Mol. Opt. 35, 5105 (2002). 

[16] A. Newell and J. Moloney, Nonlinear Optics (Addison- 
Wesley, Redwood City, CA, 1992). 

[17] T. Dohnal, D. Pelinovsky, and G. Schneider, "Coupled- 
mode equations and gap solitons in a two-dimensional 
nonlinear elliptic problem with a separable periodic po- 
tential" (to be submitted). 

[18] W. V. Houston, Phys. Rev. 57, 184 (1940). 

[19] L. D. Landau, Phys. Z. Sowjetunion 2, 46 (1932); C. 
Zener, Proc. R. Soc. London A 137, 696 (1932); E. Ma- 
jorana, Nuovo Cimento 9, 43 (1932). 

[20] S. Brundobler and V. Elser, J. Phys. A: Math. Gen. 
26, 1211 (1993); A.V. Shytov, Phys. Rev. A 70, 052708 

(2004) . 



